﻿#pragma once
#include <vector>
#include <eigen3/Eigen/Dense>
#include "Function.h"
using namespace std;
using namespace Eigen;

template <int order>
class PPspline;

// 分段线性样条插值
template<>
class PPspline<1> {
private:
    double start;                 // 插值范围起始点
    double end;                   // 插值范围结束点
    int N;                        // 插值点数量

    vector<double> values;   // 插值点对应的函数值
    int condition;

public:
    // 公有构造函数，用于创建 PPspline 的实例
    PPspline(Function& _A, double _start, double _end, int _N, int _condition) : start(_start), end(_end), N(_N), condition(_condition) {
        
        values.resize(N);
        // 计算插值点对应的函数值
        for (int i = 0; i < N; i++) {
            // 确保能够合法访问 Function 对象的 operator() 函数
            values[i] = _A(start + i * (end - start) / (N - 1));
        }
        for (int i = 0; i < N - 1; i++)
        {
            int x1 = start + i * (end - start) / (N - 1);
            int x2 = start + (i + 1) * (end - start) / (N - 1);
            int y1 = values[i];
            int y2 = values[i + 1];

            vector<double> temp;
            temp.push_back(y1 - (y2 - y1) / (x2 - x1) * x1);
            temp.push_back((y2 - y1) / (x2 - x1));
        }
    }

};

template<>
class PPspline<3> {
private:
    double start;                 // 插值范围起始点
    double end;                   // 插值范围结束点
    int N;                        // 插值点数量
    vector<double> values;   // 插值点对应的函数值
    vector<double> knots; 
    int condition;
    int n;
    vector<double> result;
    vector<double> xInterval;

public:
    PPspline(Function& _A, double _start, double _end, int _N, int _condition) : start(_start), end(_end), N(_N), condition(_condition) {
        n = N - 1;
        knots.resize(N);
        values.resize(N);
        for (int i = 0; i < N; i++) {
            knots[i] = start + i * (end - start) / (N - 1);
            values[i] = _A(knots[i]);
        }

        for (int i = 0; i < n; i++)
        {
            xInterval.push_back(knots[i + 1] - knots[i]);
        }
        if (condition == 0) {
            double m_a = _A.diff(knots[0]);
            double m_b = _A.diff(knots[N - 1]);
            complete(m_a, m_b);
        }
        else if (condition == 1) {
            double m_a = _A.diff2(knots[0]);
            double m_b = _A.diff2(knots[N - 1]);
            D2(m_a, m_b);
        }
        else if (condition == 2) {
            double m_a = 0;
            double m_b = 0;
            D2(m_a, m_b);
        }

    }

    void complete(double m_a, double m_b) {
        MatrixXd coeMatrix = calculateMatrixA_com();
        VectorXd rhs = calculateVectorC_com(m_a, m_b);
        VectorXd res = coeMatrix.lu().solve(rhs);
        for (int i = 0; i < n + 1; i++) {
            result.push_back(res(i));
        }
    }

    void D2(double m_0, double m_n) {
        MatrixXd coeMatrix = MatrixXd::Zero(n - 1, n - 1);
        vector<double> miu, lamda;
        VectorXd rhs(n - 1);

        for (int i = 1; i < n; i++) {
            miu.push_back(xInterval[i - 1] / (xInterval[i - 1] + xInterval[i]));
            lamda.push_back(1 - miu[i - 1]);
        }
        
        for (int i = 1; i < n; i++) {
            rhs(i - 1) = (6.0 * ((values[i + 1] - values[i]) / xInterval[i] - (values[i] - values[i - 1]) / xInterval[i - 1]) / (knots[i + 1] - knots[i - 1]));
        }

        rhs(0) -= miu[0] * m_0;
        rhs(n - 2) -= lamda[n - 2] * m_n;

        for (int i = 0; i < n - 1; i++) {
            coeMatrix(i, i) = 2;
        }

        for (int i = 0; i < n - 2; i++) {
            coeMatrix(i, i + 1) = lamda[i];
            coeMatrix(i + 1, i) = miu[i + 1];
        }

        VectorXd res = coeMatrix.lu().solve(rhs);

        result.push_back(m_0);
        for (int i = 0; i < n - 1; i++) {
            result.push_back(res(i));
        }
        result.push_back(m_n);
    }

    double get_value(double x) {
        int intervalIndex = findIntervalIndex(x);

        double h_temp = xInterval[intervalIndex];
        double x_minus_knots_temp = x - knots[intervalIndex];
        double x_minus_knots_temp1 = x - knots[intervalIndex + 1];

        double term1 = result[intervalIndex + 1] * pow(x_minus_knots_temp, 3) / (6.0 * h_temp);
        double term2 = result[intervalIndex] * pow(x_minus_knots_temp1, 3) / (6.0 * h_temp);
        double term3 = (values[intervalIndex] - result[intervalIndex] * h_temp * h_temp / 6.0) * x_minus_knots_temp1 / h_temp;
        double term4 = (values[intervalIndex + 1] - result[intervalIndex + 1] * h_temp * h_temp / 6.0) * x_minus_knots_temp / h_temp;

        double s = term1 - term2 - term3 + term4;
        return s;
    }

    int findIntervalIndex(double x) {
        for (int i = 0; i < n; i++) {
            if (knots[i] <= x && knots[i + 1] >= x) {
                return i;
            }
        }
        return -1;
    }

    MatrixXd calculateMatrixA_com() {
        MatrixXd coeMatrix = MatrixXd::Zero(n + 1, n + 1);

        coeMatrix.diagonal() = VectorXd::Constant(n + 1, 2);

        vector<double> miu, lamda;
        for (int i = 1; i < n; i++) {
            miu.push_back(xInterval[i - 1] / (xInterval[i - 1] + xInterval[i]));
            lamda.push_back(1 - miu[i - 1]);
        }

        for (int i = 0; i < n - 1; i++) {
            coeMatrix(i + 1, i + 2) = lamda[i];
            coeMatrix(i + 1, i) = miu[i];
        }
        coeMatrix(0, 1) = 1;
        coeMatrix(n, n - 1) = 1;

        return coeMatrix;
    }

    VectorXd calculateVectorC_com(double m_a, double m_b) {
        VectorXd rhs = VectorXd(n + 1);

        rhs(0) = (6.0 / xInterval[0] * ((values[1] - values[0]) / xInterval[0] - m_a));
        for (int i = 1; i < n; i++) {
            rhs(i) = (6.0 * ((values[i + 1] - values[i]) / xInterval[i] - (values[i] - values[i - 1]) / xInterval[i - 1]) / (knots[i + 1] - knots[i - 1]));
        }
        rhs(n) = (6.0 / xInterval[n - 1] * (m_b - (values[n] - values[n - 1]) / xInterval[n - 1]));

        return rhs;
    }
    
};